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Section  1 
INTRODUCTION 

In  July  1986,  a  research  program  was  initiated  by  S-CU3ED  and  its 
subcontractor,  ENDOCHRONICS,  INC.,  for  the  purpose  of  developing  an 
advanced  constitutive  model  of  I8ST  soils  based  upon  the  concepts  of 
endochronic  plasticity.'1*3'  The  research  was  sponsored  by  the  U.S.  Army 
Engineer,  Waterways  Experiment  Station  (WES)  as  part  of  its  continuing  effort  to 
pursue  the  development  and  validation  of  new  and  improved  analytical  methods 
for  predicting  explosively-induced  ground  motion.  In  order  to  have  a  complete 
and  reproducible  set  of  ISST  data  appropriate  for  this  purpose,  WES  performed 
a  series  of  laboratory  tests  on  reconstituted  ISST  soil  samples. w  Some  of  the 
tests  in  this  series  were  of  a  special  nature,  such  as  shearing  at  constant 
hydrostatic  pressure,  to  provide  the  type  of  data  that  are  most  convenient  for 
fitting  the  endochronic  model. 

The  model  development  was  conducted  in  two  phases,  for  reasons  that  will 
become  apparent  in  the  sequel.  In  the  first  phase,  which  was  initiated  in  July 
1986,  the  then  existing  version  of  the  endochronic  model  (see  Ref.  5,  for 
example)  was  applied  to  the  ISST  data  provided  by  WES.  It  was  realized  at  the 
time  that  this  endochronic  mode!  was  only  capable  of  predicting  compaction 
during  shearing  at  constant  hydrostatic  pressure.  Thus,  dilatancy  cannot  be 
described  by  this  model.  A  cursory  review  of  data  from  other  sandy  soils  similar 
to  the  ISST  soil  conducted  prior  to  the  program  had  indicated  that  this  limitation 
of  the  model  would  probably  not  be  serious.  The  issue  turned  out  to  be  much 
more  complex  than  we  expected,  however,  for  reasons  to  be  discussed  in  the 
text. 


Prior  to  applying  the  model  to  the  WES  data,  the  data  were  processed  on 
the  basis  of  the  usual  recommendations  of  WES,  as  described  in  Ref.  6,  that  is, 
volumetric  strains  for  cases  of  pure  hydrostatic  compression  were  defined  on 
the  basis  of  the  "uniform"  approximation,  while  volumetric  strains  generated 
during  a  process  where  shear  was  involved  were  defined  according  to  the 
"cone"  approximation.  It  was  found  that,  when  data  involving  shear  are 
processed  on  the  basis  of  the  "uniform”  approximation,  dilatancy  is  usually 
observed  while,  if  the  "cone"  approximation  is  used  to  process  the  data, 
compaction  will  be  predicted.  Thus,  there  is  the  potential  for  serious 
inconsistencies  between  predictions  and  data,  which  were  vividly  confirmed  in 
the  first  phase  of  this  study. 
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The  inconsistency  arises  due  to  the  fact  that  the  purely  hydrostatic 
component  of  a  model  is  usually  fit  to  data  which  has  been  processed  on  the 
basis  of  the  uniform  approximation.  This  component  is  then  used  as  part  of  the 
more  general  model  to  predict  the  volumetric  component  of  response  in  cases 
of  non*isotropic  loading,  involving  shear.  In  such  cases,  the  predicted 
volumetric  response  will  be  the  result  of  a  model  based  on  the  uniform 
approximation,  while  the  data,  processed  in  the  recommended  fashion,  would 
be  based  on  the  "cone”  approximation.  Because  of  this,  the  model  will  most 
likely  predict  dilatancy  while  the  data  will  probably  show  compaction.  The 
results  obtained  from  the  first  phase  of  this  study  in  which  the  endochronic 
model  was  fit  to  the  ISST  data  in  the  recommended  manner  dramatically 
confirmed  these  inconsistencies. 

Unfortunately,  there  appears  to  be  no  precise  way  of  accurately  defining 
the  volumetric  strain,  needed  for  constitutive  model  development,  from  the  usual 
triaxial  data.  The  difficulty  arises  from  the  fact  that  the  strain  fields  in  triaxial  tests 
ore  typically  inhomogeneous.  In  view  of  this,  and  given  that  the  usual  data  from 
the  triaxial  test  gives  the  nominal  engineering  axial  strain  e  over  the  specimen 
length  and  the  nominal  engineering  radial  strain  er  at  the  midpoint  of  the 
specimen,  it  appears  that,  while  admittedly  nor  rigorous,  the  simple  expression 


e  *  eg  2er  (1) 

appears  to  provide  the  best  definition  of  the  nominal  engineering  volumetric 
strain  that  is  possible  within  the  limitations  of  the  data  provided.  Several 
techniques  for  more  accurately  measuring  strain  in  triaxial  tests  have  recently 
been  proposed  (see  Refs.  7  and  8,  for  example).  Until  such  techniques  have 
been  implemented;  however,  Eq.  (1)  appears  to  be  the  preferred  way  of  defining 
volumetric  strain  from  triaxial  data. 

In  view  of  the  difficulties  encountered  in  the  first  phase  of  the  study  due  to 
the  use  of  the  uniform  and  cone  approximations,  it  was  mutually  agreed  by  WES 
and  S-CUBED  that  a  second  phase  should  be  undertaken,  using  data  that  are 
processed  only  on  the  basis  of  Eq.  (1),  which  to  first  order  is  the  uniform 
approximatior  This  eliminates  the  inconsistencies  that  can  arise  from  the  joint 
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use  of  the  "cone"  and  "uniform"  approximations.  However,  another  problem 
develops,  since  the  use  of  Eq.  (1)  in  conjunction  with  the  ISST  data  will  lead  to 
dilatancy,  which  is  not  within  the  scope  of  the  endochronic  model  used  in  the 
first  phase.  Therefore,  an  endochronic  model  with  the  capability  to  describe 
dilatancyis  needed.  Fortunately,  such  a  model  had  very  recently  become 
available, but  had  not  been  val, dated  or  applied  to  real  materials. 

The  second  phase  of  this  study,  which  began  in  April  1987,  was 
undertaken  for  the^purpose  of  exploring  the  application  of  the  new  dilatant 
endochronic  model™’10'  to  the  ISST  soil  data,  using  Eq.  (1)  to  define  volumetric 
strain  from  the  data.  The  present  report  has  been  prepared  mainly  to 
summarize  this  work  In  Section  2,  the  new  dilatant  endochronic  model  is 
described  and  the  special  case  of  shear  at  constant  pressure  is  considered.  The 
application  of  the  model  to  the  ISST  data  is  described  in  Section  3.  A  numerical 
scheme  and  corresponding  computer  program  for  computationally  dealing  with 
the  model  are  described  in  Section  4.  Also,  in  this  section,  the  case  of  shear  at 
constant  pressure  is  considered  numerically.  Section  5  discusses  some 
difficulties  that  were  encountered  in  the  numerical  study,  and  suggests  possible 
ways  to  solve  them.  Finally,  the  conclusions  drawn  from  this  study,  as  well  as 
recommendations  for  further  study,  are  given  in  Section  6. 
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AN  ENDOCHRONIC  SOIL  MODEL  WITH  DILATANCY 
2.1  BASIC  CONSIDERATIONS. 

The  application  of  endochronic  plasticity  to  soils  has.  so  far,  been  limited 
either  to  one-dimensional  histories  (shear  and  hydrostatic  response)  under 
cyclic  and  other  loading-unloading  conditions  or  to  two-dimensional  histories 
(shear-volume  interaction)  under  densifying  conditions.  The  question  of  dilatant 
behavior  has  hitherto  remained  unresolved,  in  the  sense  that,  until  now,  the 
constitutive  relation  in  question  [11]  gave  rise  to  densification  in  the  presence  of 
shearing  under  constant  hydrostatic  stress. 

What  is  needed  is  a  constitutive  equation  which  under  conditions  of  low 
density  and/or  high  pressure  will  give  rise  to  shear-induced  densification  but 
under  conditions  of  high  density  and/or  low  pressure  will  give  rise  to  shear- 
induced  dilatancy.  In  this  section  we  present  a  thermodynamic  approach  which 
leads  to  such  characteristics  of  soil  behavior.  The  full  analysis  is  given  in  detail 
in  Ref.  10. 

The  reasoning  that  led  to  the  present  treatment  is  broadly  as  follows.  The 
coupling  between  deviatoric  and  hydrostatic  behavior,  that  ultimately  leads  to 
dilatant  deformation,  must  come  from  three  possible  sources: 

(i)  The  intrinsic  time  through  the  shear-volumetric  coupling 
parameter  k; 

(ii)  The  expression  for  the  free  energy; 

(iii)  The  rate  equations  for  the  internal  variables. 

Source  (i)  alone,  will  always  give  rise  to  densification,  as  the  application  of  the 
relevant  equations  actually  showed.  Source  (ii)  is  not  physically  viable  because 
given  a  soil  with  a  certain  porosity  the  onset  of  dilatancy  under  monotonic 
shearing  is  governed  by  the  prevailing  hydrostatic  stress.  If  (ii)  is  to  be  the 
source,  then  a  change  in  the  form  of  the  free  energy  must  take  place  upon 
varying  the  hydrostatic  stress,  a  phenomenon  which  does  not  appear  physically 
plausible.  One  would  expect  the  form  of  the  free  energy  to  remain  invariant  with 
a  change  in  the  hydrostatic  stress.  The  remaining  plausible  cause  is  source  (iii) 
and  this  is  the  one  that  we  developed  in  Ref.  10  and  summarize  in  the  present 
work. 
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2.2  FORMULATION  OF  MODEL 

In  Ref.  10  took  the  position  that  whereas  the  deviatoric  plastic  work  is  a 
cause  of  hydrostatic  plastic  strain,  it  is  external  to  the  hydrostatic  process  in  the 
sense  that  it  is  not  a  hydrostatic  mechanism.  Therefore,  it  qualifies  as  a 
thermodynamic  internal  force  of  the  first  kind  in  the  sense  of  Ref.  12.  In  the 
presence  of  such  forces  the  thermodynamic  equations  appropriate  to  the  rigid 
plastic  solid  that  represents  the  plastic  behavior  of  the  soil  are  the  following,  in 
the  usual  notation,  in  terms  the  deviatoric  and  hydrostatic  free  energies  tip  and 
,'l_l  respectively.  For  a  more  detailed  thermodynamic  treatment  see  Ref.  13. 


to  =  5  y~ A:'  I  |sP  ■  3r||2 
to  =  \  y  Br  l‘p  ■  qr|2 


8to  ,  r  dSr 

Wr  +  an  ar 


=  0 


r  *>r 

22  ar 


where 


•  D 

dz 


(2) 

(3) 

(4a,  b) 

(5) 

(6) 

(7) 
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Thus  Rf  is  a  thermodynamic  force  of  the  first  kind,  i.e.,  one  which  is 
internally  applied  but  which  is  external  to  the  (internal)  mechanism  on  which  it 
operates.  The  role  of  the  coefficient  is  to  determine  the  degree  to  which 
the  plastic  work  affects  the  r’th  hydrostatic  mechanism. 

The  intrinsic  time  dz  is  given  by  the  customary  relation 


dZ2»||d8P||2  +  k2ldeP|2  (8) 

The  resistance  coefficients  eL  1  and  alg  are  not  constant  but  are  related  to  the 
hardening  functions  Fq  and  F|_|  respectively  by  the  equations 


—  Fr»  a, 


ro 


D  “11 


where 

FD=  FD<ff  Z)  '  FH=FH(eP) 


(9a,  b) 


(10a, b) 


The  dependence  of  Fq  on  z  is  weak  so  that  for  deformations  other  than  cyclic, 
for  which  the  variation  in  z  is  small,  the  latter  plays  a  minor  role  in  Fq  and  may 
be  ignored.  The  question  as  to  whether  the  dilatancv  coefficients  ag-j  are 
constant  during  the  deformation,  or  whether  they  change  with  hardening,  is  an 
open  one  at  the  moment.  In  this  work  they  have  been  taken  as  constant  and  this 
choice  seems  to  be  consistent  with  experiment. 

A  straightforward  analysis  using  Eqs.  (2)  to  (7)  in  light  of  the  initial 
conditions 


Qr(0)-Q.  qr(0)=0 


(ila, b) 


gives  rise  to  the  following  set  of  two  constitutive  equations  that  govern, 
respectively,  the  deviatoric  and  hydrostatic  behavior  of  soils: 


6 


SSS-R-88-9375 


•  z')ar dz' 


(12) 


a  - 


'Ph  •  z')ar  *'  * 


dgp 

f(2h  -  r)  §  •  ard2' 


(13) 


*'0  *'0 

where  the  kernels  p,  +  and  r  are  all  weakly  singular  but  integrable  in  any  finite 
domain  [0,z]  and: 


_j_  _  dz  .  _  dz 


(14a, b) 


P  = 


Ar  e 


(15) 


*•  21  Br 


B,  e 


-*r*H 


(16) 


a21  e’#r1H 


(17) 


=  Aj_  =  Bj_ 

*r  „ro  ’  ^r  “  „ro 


a 
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®22 


(18a, b) 


Some  simplifications  occur  when  the  effect  of  the  deviatoric  plastic  work  rate  is 
distributed  uniformly  among  the  hydrostatic  mechanisms,  in  which  case 
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a21  =  a21  (for  al  1  r)  (19) 

In  this  event,  and  in  view  of  Equations  (16)  and  (17) 


rM=a2l'(zH)  ™ 

so  that  once  ^(zH)  is  known,  from  a  simple  monotonic  hydrostatic  test,  then 
r(zH)  is  also  determined  within  a  multiplicative  constant. 

In  addition  there  exist  the  elastic  relations 


dge  =  d§/2G  ;  de®  =  da  IK  (21) 

where  /i  is  the  elastic  shear  modulus,  and  K  the  elastic  bulk  modulus  neither  of 
which  need  be  constant  in  the  course  of  deformation.  The  total  strain  is  given  by 
the  relations: 


dg  =  dg®  +  dgp  ;  de  *  de®  +  dep 


(22) 


The  general  method  for  determining  the  form  of  the  functions  p,  <p  and  r,  the 
functions  FH  and  Fq,  the  elastic  moduli  G  and  K  as  well  as  the  parameter  k, 
from  appropriate  experiments,  is  a  subject  that  will  be  discussed  later  on  in  this 
section.  The  specific  determination  of  these  functions  and  parameters  for  the 
WES  ISST  soil  is  addressed  in  Section  3. 

2.3  GENERAL  APPROACH  FOR  DETERMINING  MATERIAL  FUNCTIONS 

AND  PARAMETERS. 

In  this  section,  the  general  approach  for  determining  the  material  functions 
and  parameters  of  the  endochronic  model  described  in  the  preceding  section  is 
given.  The  functions  and  parameters  that  need  to  be  determined  for  a  particular 
material  are:  K,  G,  k,  ^(zD),  ^(z^)*  fh  and  FD’ 
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The  bulk  modulus  K,  hydrostatic  kernel  #(z^)  and  the  hardening  function 
Fh  can  be  determined  from  a  pure  hydrostatic  test,  which  should  have  at  least 
two  unloadings  and  reloadings  and  extend  well  into  the  concave  part  of  the 
virgin  hydrostatic  stress-strain  curve.  The  data  should  consist  of  continuous 
measurements  of  a  and  e. 

The  shear  modulus  G,  coupling  parameter  k,  dilatancy  parameter  rQ, 
hardening  function  Fq  and  the  shear  kernel  p(zD)  are  most  efficiently  determined 
from  triaxial  tests  in  which  specimens  are  sheared  under  fixed  hydrostatic  (not 
confining)  pressure.  The  fixed  hydrostatic  pressures  selected  for  these  tests 
should  adequately  cover  the  hydrostatic  range  of  interest;  they  should  be 
reached  monotonically  and  lie  on  the  concave  part  of  the  virgin  hydrostatic 
curve.  The  shearing  should  be  taken  out  to  where  the  shear  stress  has 
essentially  reached  a  limiting  value.  The  data  should  consist  of  continuous 
measurement  of  o,  e,  §  and  §. 

Specific  details  of  the  approach  for  determining  the  material  functions  and 
parameters  from  the  types  of  data  discussed  above  are  given  below. 

2.3.1  The  Bulk  Modulus  K. 

The  bulk  modulus  K  at  the  onset  of  hydrostatic  deformation  is  determined 
by  the  slope  of  the  hydrostatic  stress-strain  curve  at  the  origin.  We  denote  this 
value  of  K  by  KQ.  However,  in  the  case  of  ISST  soils  it  was  found  that  K  does  not 
remain  constant  but  varies  with  compaction.  Simultaneously  it  becomes 
dependent  on  the  elastic  hydrostatic  strain.  This  complex  behavior  is  inferred  by 
making  the  following  observations: 

(a)  Upon  loading  to  some  value  of  a  and  unloading  to  zero  stress  and 
reloading  one  finds  that  K  (i.e.,  the  slope  of  the  hydrostatic  stress-strain  curve  at 
the  onset  of  reloading)  is  not  equal  to  KQ. 

(b)  The  reloading  stress-strain  curve  almost  retraces  the  unloading 
curve  with  the  implication  that  unloading  is  virtually  elastic,  and  because  it  is 
non-linear,  obviously  dependent  on  the  elastic  strain.  This  also  implies  that 
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plastic  unloading  in  the  stress-plastic  strain  space  is  almost  vertical  with  the 
inference  that  ^(zH)  is  close  to  a  delta  function.  In  mathematical  terms  K  admits 
the  representation: 


K~K(ee,€P)  (23) 

In  ISST  soils  Eq.  (23)  is  closely  approximated  by  the  multiplicative  form 
given  in  Eq.  (24),  i.e., 


K-KjePiKjU®)  (24) 

where,  without  loss  of  generality,  we  may  set  K1  (0)  =  1. 

Now  let  P  be  a  point  on  the  hydrostatic  strain  axis,  reached  upon  unloading  from 
a  value  of  a.  Then  KQ  is  determined  from  the  initial  slope  of  the  reloading  hydrostatic 
stress-strain  curve  at  P  (where  e®  =  0)  and  K1  from  the  shape  of  the  unloading  part  of 
the  hydrostatic  curve,  terminating  at  P. 

2.3.2  The  Hardening  Function  F^. 

At  this  point  certain  observations  are  in  order.  The  hydrostatic  response  under 
purely  hydrostatic  conditions  is  given  by  Eq.  (25),  obtained  by  setting  §  equal  to  zero  in 
Eq.  (13).  Thus 
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If  Fu  were  constant  and  equal  to  unity,  say,  then  in  view  of  Eqs.  (25)  and  (26),  z». 
would  be  equal  to  «P  and  Eq.  (25)  would  give  a  as  a  convex  function  of  e*5. 
However  in  actual  experiments  a  becomes  asymptotically  a  very  rapidly 
increasing  concave  function  of  ep  as  ep  increases.  In  particular  this  function  is 
invariably  of  an  exponential  type,  so  that 


o  -  #0  expUep) 


(27) 


for  large  values  of  ep.  Equation  (27)  is  a  basic  hypothesis  of  the  critical  state 
theory.  If  were  a  delta  function  in  the  sense  of  Eq.  (28),  then, 


* W  88  *0  'PU 


(28) 


Thus  in  view  of  Eqs.  (25),  (8)  and  (14b): 


a 


(29) 


Comparison  of  Eqs.  (27)  and  (29)  yields  the  important  result, 


F|_j  =  exp(/Jcp)  (30) 

Consider  now  the  other  simple  case  where  the  hydrostatic  response  is  given  in 
terms  of  one  internal  variable.  In  this  event  the  resulting  constitutive  response  is 
given  by  the  following  equation: 


A  +  (da/dz^  dcP/dZn) 


(31) 


which  under  monotonic  loading  conditions  becomes  the  nonlinear  differential 
equation: 
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A  +  (d»/d«P)FH  -  rt0  Fh  ,  (32) 

where  p*  and  pQ  are  material  constants. 

We  now  ask  the  question  as  to  what  should  the  form  of  the  hardening 
function  be  if  the  hydrostatic  response  is  to  be  given,  albeit  asymptotically  by 
Eq.  (27).  A  simple  substitution  for  o  (as  given  by  Eq.  (27))  in  Eq.  (31)  gives  rise 
to  the  following  result: 


Fh  =  exp(^ep)/{l  -  a  exp(^ep)}  (33) 

where  a  -  pip*.  Now  FH  is  more  complicated  but  still  looks  a  great  deal  like  an 
exponential  function  for  small  values  of  "a".  It  is  of  interest  to  note  that  in  view  of 
Eq.  (32)  it  follows  that  as  p *  *  •,  ^(zJ  tends  to  a  delta  function,  a  is  given  by 
Eq.  (29)  and  FH  by  Eq.  (30),  which  is  the  limit  of  Eq.  (33)  as  a  +  0. 

2.3.3  The  Kernel  Function  d(z^). 

The  function  f  (zH)  can  be  determined  uniquely  from  Eq.  (25)  if  the  form  of 
the  hardening  function  Fu  is  known  a  priori.  To  that  end  we  note  that  there  is 
sufficient  experimental  evidence  to  support  the  position  the  FH  is  of  exponential 
character  and  is  given  by  Eq.  (30).  Under  conditions  of  monotonic  loading 
Eq.  (25)  then  becomes: 


a 


f*H 

*Ph  '  z')FH<z'>dz' 
Jo 


(34) 


where 


"W  :  fhW  ~fh{cP(zh)} 


(35a, b) 
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in  the  case  where  Eq.  (30)  applies  a  simple  calculation  shows  that 


fh"1/(1*^!h)  (36) 

Since  now  both  a  and  FH  are  known.  Eq.  (34),  which  is  a  Volterra  integral 
equation,  may  be  solved  numerically  to  give  the  unknown  function  f(zH). 

In  the  actual  case  of  ISST  soil,  the  steepness  of  the  unloading  hydrostatic 
stress-strain  curve  indicated  that  l(z^)  is  indeed  dose  to  a  delta  function.  In  this 
event  the  numerical  values  of  l  are  so  dominant  near  zH  *  z'  that,  by  virtue  of 
the  mean  value  theorem,  one  can  represent  the  integral  on  the  right-hand  side  of 
Eq.  (34)  by  the  expression: 


a 


-  2-)dz- 


'o 


A 


rH 


#(z*)ck* 


(37) 


Thus,  in  view  of  Eqs.  (35)  and  (36): 


-  pwwf'ftip  -  Phi)  m 

2.3.4  The  Material  Constants  k  and  rQ. 

We  begin  with  the  hydrostatic  constitutive  equation,  previously  Eq.  (13), 
which  in  the  case  of  the  cylindrical  triaxial  test  becomes: 


A 


a  = 


SH 


'fa  ’  Z')$F  dz'  + 


rfa  ’  z')  s  *  Hr  dz' 


(39) 
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where  s  -  J*5 75(*2  •  a..)  end  e^  «  f£75(«§  *  «?),  with  x1  being  the  exial 
direction. 

The  sheer  modulus  Q  is  presumed  known,  its  determination  will  be 
discussed  later  in  the  section.  Thus  e^  may  be  calculated  from  e  via  Eqs.  (21) 
and  (22). 

Previously  f  and  r  were  represented  by  delta  functions  in  the  specific  case 
where  initial  hydrostatic  stress  lay  on  the  concave  part  of  the  hydrostat.  Here  we 
consider  the  more  general  case  where  p  and  T  are  represented  by  an 
exponential  function.  Thus  we  set 


4 


T  =  p*TQB 


(40) 


Note  that,  in  the  limit  as  the  rignt-hand  sides  of  Eqs.  (40a, b)  become  delta 
functions,  thereby  recovering  the  previous  case. 

Using  Laplace  transforms  one  may  now  solve  Eq.  (39)  in  the  light  of 
Eqs.  (41)  and  obtain  the  following  expression  for  e*5: 


(41) 


Now  let  a  ~  a  ,  where  o  is  the  constant  hydrostatic  stress  at  which 
shearing  occurs,  ana  let  dc£  be  the  increment  in  volumetric  strain  due  to  shear, 
at  constant  hydrostat'C  stress.  Then  differentiating  Eq.  (41)  we  obtain  the 
following  relation  for  de^: 


HI 

=  ao  *  ro 

Is  •  arJ 

(42) 
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Note  that  the  same  equations  would  have  been  obtained  had  we 
represented  f  and  r  by  ^-functions.  However,  now  there  is  an  important 
difference  in  that  *Q  is  no  longer  equal  to  e^€P  where  «p  is  the  volumetric  strain 
at  the  end  of  the  (initial)  hydrostatic  phase.  In  met,  in  the  present  case: 


f*. 

» -  p  p0»  («3) 

‘'o 

or 

^  (**) 

The  effect  is  illustrated  in  Figure  1 .  Note  that  at  low  stresses  the  form 
l  -  *0ff(zH)  seriously  overestimates  the  initial  hydrostatic  stress. 

let  us  now  recall  Eqs.  (8)  and  (42)  with  the  constraint  a  «  aQ,  i.e., 


dz2  =  |  \6sP\  |2  +  ^(deP)2  (45) 

=  *0e^ekdep  +  •  dgp  (46) 

except  that  now  aQ  ?  e^c  but  is  given  by  Eq.  (43).  Specifically  in  the  case  of  the 
triaxial  test: 

dz2  =  (dep)2  +  k2tdep)2  (47) 

oQdz  =  ^0e^6  kdep  +  rQs  dep  (48) 
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W«r*  *-{§(«,  •  *2)  »nd  tP  -  J1  («?  •  «|)  (48-) 

Equations  (47)  and  (49)  can  be  solved  simultaneously  to  give: 


The  constant  r  is  determined  at  the  point  A  where  the  slope  of  the  curve 
of  the  shear-induceanydroststic  strain  versus  plastic  shear  strain  is  zero.  (See 
Fig.  2).  Then  it  follows  from  Eq.  (49)  that 


r 


o 


(50) 


where  s»  is  the  value  of  s  at  point  A  in  Fig.  2.  Since  the  right-hand  side  of 
Eq.  (49)  is  now  known,  defyde*5  may  be  found  for  various  values  of  k.  The  one 
which  provides  the  best  agreement  with  the  experimental  data  is  chosen. 

Approximate  Solution  of  Equation  (49). 

The  hydrostatic  response  can  be  determined  explicitly  by  means  of  an 
approximate  solution  of  Eq.  (43).  In  effect,  since  p*  exp(-0*zH)  is  close  to  a  delta 
function,  Eq.  (43)  becomes: 


* 


Thus 


(51) 
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(52) 


Equation  (52)  may  now  be  used  in  Eq.  (50)  to  give  the  following  explicit  result 


and 


(55) 


Solution  of  Eq.  (49)  for  Asymptotically  Small  aQ. 

Of  interest  is  the  case  of  asymptotically  small  aQ  and  e£  in  which  event 


- 


(56) 


and 
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1 


(57) 


Thus  Eq.  (54)  becomes 


D  = 


(58) 


It  follows,  therefore,  that  at  small  values  of  aQ,  D  is  magnified  and  thus  the  value 
of  dcR/dep  is  depressed  and  the  onset  of  dnatancy  is  accelerated  in  accordance 
with  ooservation. 

2.3.5  The  Shear  Modulus  G. 

The  singularity  of  the  shear  kernel  ensures  that  the  behavior  is  always 
elastic  at  loading,  unloading  and  reloading  points,  on  the  shear  stress-strain 
curve.  Thus  the  initial  shear  modulus  is  determined  by  the  slope  of  the  tangent 
at  the  origin  of  the  shear  stress-strain  curve.  Measurement  of  this  slope  at  other 
loading  and  unloading  points  will  determine  if  G  is  constant.  If  not,  then  G  will 
most  likely  depend  on  the  second  invariant  of  the  elastic  deviatoric  shear  strain 
tensor,  i.e., 


G  =  G 1 1  IS 


.e 


(59) 


as  was  found  to  be  the  case  in  plain  concrete  (see  Ref.  5).  The  form  of  the 
relationship  will  be  found  by  measuring  G  at  such  points.  In  the  present  study  G 
is  considered  constant  and  was  determined  from  the  initial  slope  of  the  shear 
stress-strain  curve. 

2.3.6  The  Deviatoric  Memory  Kernel  j(zp). 

The  kernel  p  is  determined  from  the  deviatoric  response  equation  (12),  i.e., 
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S 


A 

deP 

'Pd  *  Z1  Sb  d2' 


(60) 


at  a  hydrostatic  stress  value  oQ  at  which  the  deviatoric  hardening  function  Fq  is 
normalized  and  equal  to  unity,  i.e., 


fdW  - 1 


(61) 


The  determination  of  Fq  for  other  values  of  a  will  be  discussed  shortly.  Thus  in 
the  course  of  the  sheamg  test  conducted  ala  =  oQ,  we  have 


dzjj^dz 


(62) 


Now  let  z  be  the  value  of  z  at  the  conclusion  of  the  hydrostatic  test.  Then  in 
view  of  Eq.  (62) 


h>  “ z  •  zo  I63) 

Also  let  the  right-hand  side  of  Eq.  (49)  be  denoted  by  a  function  7(6^).  This 
function  is  known  from  experimental  measurement.  Thus 


Combining  Eqs.  (62)  and  (47)  one  finds  the  relation 


dep 

®B  17 


77P75 


(64) 
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Integration  of  Eq.  (65)  yields  the  relation 

ep  =  ePfo) 

Let 

dep  dfif  _r_  x 

"  Opbi) 

where  g  is  a  known  function  of  Zp.  Also 

s(ep)  =  s[eP(zD]]  -  sfzj 
In  view  of  Eq.  (60) 


(66) 


(67) 


(68) 


sPo) 


P[z d  -  z‘)g(z')dz' 

u 


(69) 


Equation  (69)  is  a  Volterra  integral  equation  of  the  first  kind  which  can  be  solved 
for  p(zD)  since  both  s  and  g  are  known.  The  method  of  solution  is  given  in 
Ref.  3. 

2.3.7  Determination  of  the  Hardening  Function  Fp. 

The  analysis  behind  the  procedure  v/as  discussed  at  length  in  Ref.  5.  Let 
s„  be  the  value  of  s  in  the  limit  of  large  ep,  in  the  course  of  a  shear  test  at 
constant  hydrostatic  stress  aQ.  Then  it  was  shown  in  Ref.  5  that 


fd('o)  -  r0  <70> 

As  a  result,  the  dependence  of  FD  on  aQ,  can  be  determined  from  data  on 
shear  tests  conducted  at  several  values  of  fixed  hydrostatic  pressure  aQ  which 
span  the  hydrostatic  range  of  interest. 
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Section  3 

APPLICATION  TO  ISST  SOIL  DATA 

In  this  section,  the  endochronic  model  described  in  Section  2.2  is  applied 
to  ISST  soil,  using  the  laboratory  data  for  reconstitutied  ISST  soil  generated  by 
WES  specifically  for  this  study  [4].  Before  the  model  w^s  fitted,  the  WES  data 
were  transformed  from  engineering  strain  to  natural  strain  so  that  the  resulting 
model  would  be  appropriate  for  use  in  conjunction  with  existing  wave 
propagation  codes,  which  typically  are  formulated  in  terms  of  natural  strain.  In 
addition,  the  volumetric  strain  was  determined  from  the  data  in  a  consistent 
manner,  using  Eq.  (1)  given  earlier.  The  material  functions  and  parameters  of 
the  model  were  determined  from  the  WES  ISST  data  according  to  the  methods 
described  in  Section  2.3.  Specific  details  of  the  manner  in  which  this  was 
accomplished,  as  well  as  the  results,  are  given  below. 

3.1  DETERMINATION  OF  HYDROSTATIC  FUNCTIONS  AND 
PARAMETERS. 

The  measured  behavior  of  ISST  soil  to  pure  hydrostatic  loading  is  shown  [n 
Figure  3,  where  the  results  from  three  separate  tests  have  been  plotted.  As  this 
figure  reveals,  the  unloading-reloading  paths  practically  coincide,  indicating  that 
the  behavior  along  these  paths  is  essentially  elastic.  In  view  of  this,  the 
dependence  of  the  bulx  modulus  K  on  elastic  volumetric  strain  ee  (see  Eq.  (24)) 
was  obtained  from  the  unloading  curve  emanating  from  57  MPa.  The  expression 


K  =  Kq  +  K,  (ee)m 


(71) 


where 


K^SOIVPa,  K,  «  1.426  x  1010  IVPa,  m=  5.2  (72) 

was  found  to  describe  the  data  very  well. 


* 


See  Appendix  A  for  details. 
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Figure  3. 


10  15  20 

volumetric  strain  (%) 


Isotropic  compression  of  ISST  soils,  showing  comparison 
between  endochronic  model  and  data. 
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The  hydrostatic  hardening  function  F^(ep)  was  determined  as  follows. 
First,  a  plot  of  pressure,  9,  versus  plastic  volumetric  strain,  ep,  was  constructed 
from  the  data  shown  in  Figure  3,  using  the  expression  for  the  bulk  modulus 
given  in  Eq.  (71)  to  define  the  elastic  volumetric  strain.  Since  the  unloading- 
reloading  curves  practically  coincide,  the  kernel  function  l(Zu)  must  be  very 
dose  to  a  delta-function  and,  for  the  present  purposes,  was  in  fact  assumed  to 
be  one.  When  this  is  the  case,  the  virgin  hydrostat  is  given  by  the  expression 
(see  Eq.  (29)): 


(73) 


This  was  fit  to  the  a  versus  ep  data,  using  Eq.  (30),  with  the  result 


Fh  -  e*P  (74) 

where 

p  =  19.5  (75) 

Turning  now  to  the  hydrostatic  kernel  function  #(zH),  we  adopted  the  form 
given  by  Eq.  (16),  i.e., 


N  -p  z u 

§Bre  <76> 

and  set  N  =  2.  To  determine  the  Br  and  pr  we  proceeded  as  follows.  First,  let 
us  set: 


s  'W  t1  •  *h) 


(77) 
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A  plot  was  made,  in  terms  of  9(2^)  versus  zh*  °*  the  second  reloading  path 
shown  in  Figure  3.  This  plot  shows  that  if  approaches  an  asymptotic  value,  call 
it  ijm ,  as  Zu  *  •.  A  further  plot  was  made  of  ijm  -  9(2 u)  versus  zH,  and  the 
resulting  curve  was  fit,  via  Prony’s  method  (see  Ref.  14,  for  example),  by  an 
expression  of  the  form: 


2  *4_Z u 

In  *  V(*h)  *  ^  Cre 


(78) 


Differentiating  this  with  respect  to  zH,  and  recalling  Eq.  (38),  it  follows  that 
2  -Pr*u 

=  2=  Br«  ™  (TO) 

where 

Cr  (80) 

T,;  •.  since  pT  and  C.  are  known  from  the  Prony  fit  in  Eq.  (78),  the  Br  are  also 
knc.  5  as  a  result  of  tq.  (80).  In  this  manner,  the  following  values  of  Br  and  pr 
were  -atermined: 

B1  -  1,680  MPa  ^=  2,595. 

(81) 

^  =11, 610  MPa  P2  —  13,680. 

This  completes  the  specification  of  the  hydrostatic  functions  and 
parameters  for  ISST  soil.  A  comparison  between  the  model  predictions  based 
upon  these  functions  and  parameters  and  the  corresponding  WES  soil  data  [4] 
is  shown  in  Figure  3.  As  an  inspection  of  this  figure  reveals,  the  model  does  a 
remarkably  good  job  of  capturing  all  of  the  details  of  the  data. 
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3.2  DETERMINATION  OF  SHEAR-VOLUMETRIC  COUPLINQ 

PARAMETERS  k  AND  r0. 

For  the  determination  of  k  and  r0,  we  will  use  Eq.  (49)  and  the  data  given 
in  Ref.  4  for  the  cases  of  shear  at  various  fixed  hydrostatic  pressures.  Attention 
is  specifically  focussed  on  the  three  cases  for  which  the  fixed  hydrostatic 
pressure  lies  on  the  concave  portion  of  the  virgin  hydrostat,  i.e.,  those  with  fixed 
pressure  of  3.45, 6.90  and  10.34  MPa. 

The  importance  of  Eq.  (49)  lies  in  the  fact  that  by  fitting  it  to  the  appropriate 
shear-volumetric  coupling  data  obtained  from  shear  tests  at  constant  hydrostatic 
pressure,  the  values  of  the  parameters  k  and  rQ,  can  be  found  directly.  Thus, 
the  shear-volumetric  coupling  parameters  do  not  need  to  be  determined 
iteratively.  It  will  be  recalled  that  the  parameter  fi  was  obtained  independently 
from  the  pure  hydrostatic  compression  data  (see  Eq.  (75))  for  ISST  soils  and  has 
the  value  19.5. 

It  should  be  noted  that  Eq.  (49)  was  derived  on  the  assumption  that  the 
shear  strain  increases  monotonically.  However,  the  shear  tests  at  constant 
pressure  reported  in  Ref.  4  were,  for  other  reasons,  not  performed  under 
monotonic  shearing  conditions.  In  each  test,  the  shearing  was  interrupted  at  a 
shear  strain  of  3  percent  by  an  unloading-reloading  process,  and  then  the 
shearing  was  continued  out  to  a  shear  strain  of  about  1 1  percent,  after  which 
unloading  took  place.  The  experimental  data  for  these  tests  are  shown  in 
Figures  4  to  6,  where  the  volumetric  strain  due  to  shear  has  been  plotted 
against  the  plastic  octahedral  shear  strain.  These  figures  reveal  that  the  soil  first 
undergoes  compaction  up  to  a  peak,  followed  by  dilatancy  with  further  increases 
in  the  shear  strain. 

In  applying  Eq.  (49)  to  the  data  depicted  in  Figures  4  to  6,  we  attempted  to 
compensate  for  the  effects  of  the  unloading-reloading  process  at  3  percent 
shear  strain.  A  numerical  program  was  developed  to  integrate  Eq.  (49),  and 
using  this,  the  following  values  of  k  and  rQ  were  found  to  provide  the  best 
correlation  of  the  data  from  the  three  tests: 

k  =  0.6 

(82) 

T  =  1.25  MPa. 
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Figure  6.  Shear-volumetric  coupling  at  a  fixed  hydrostatic  pressure  of 
10.34  MPa. 
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Figures  4  to  8  show  comparisons  bstwasn  the  data  and  Eq.  (49)  with  the  above 
values  of  the  parameters.  An  inspection  of  these  figures  reveals  that  the  model 
captures  die  essential  features  of  the  data  quite  well,  in  view  of  the  fact  that  the 
data  were  not  obtained  under  truly  monotonic  shearing  conditions.  These 
figures  demonstrate  the  ability  of  the  model  to  describe  both  compaction  and 
dBatancv  in  accordance  with  experimental  observation. 

DETERMINATION  OP  DEV1ATORIC  FUNCTIONS  AND  PARAMETERS. 

Figure  7  summarizes  the  measured  dependencies  of  the  octahedral  shear 
stress  on  the  octahedral  shear  strain  from  the  nine  shear  tests  reported  in  Ref.  4 
that  were  performed  at  constant  hydrostatic  pressures.  In  view  of  the 
remarkable  consistency  and  reproducibility  of  these  data,  the  care  with  which 
these  experiments  were  done  is  dearly  evident. 

The  deviatoric  functions  G,  p(Zq)  and  Fq  were  determined  from  one  of  the 
tests  depicted  in  Figure  7,  namely,  test  RDC  624  in  which  the  fixed  hyhdrostatic 
pressure  was  6.90  MPa.  The  volumetric  strain  caused  by  the  shearing  was  taken 
to  be  irreversible  (fully  plastic)  and,  in  the  developments  which  follow,  will  be 
denoted  by  efj.  From  the  data  for  test  RDC  624,  numerical  tables  of  the  relations 


r  *  rtyP)  (83) 

«f>  -  «*) 


were  constructed,  where 


f  -  7  •  5B  (85) 

Here,  r  and  7  denote,  respectively,  the  octahedral  shear  stress  and  the 
octahedral  shear  strain.  Furthermore,  is  the  plastic  component  of  7,  while  G 
is  the  shear  modulus. 
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The  shear  modulus,  and  its  dependence  on  7,  was  found  as  follows.  The 
response  of  the  soil  at  the  initiation  of  shear  loading,  as  well  as  at  the  two  points 
of  unloading  shown  in  Figure  7,  is  assumed  to  be  purely  elastic.  On  this  basis, 
the  slopes  of  the  r  versus  7  curve  at  these  points  defines  the  corresponding 
values  of  2G.  In  this  manner,  the  following  dependence  of  2G  on  7  was  found  to 
provide  a  reasonable  description  of  the  data: 


2G  =  2^,  +  corf 


where 


23^  *  500  IVPa 


cQ  =  3000  MPa 


(86) 


(87) 


To  determine  the  kernel  function  />(zD),  we  begin  by  rewriting  Eq.  (69)  in 
the  form 


r(w)  = 


f/>(w  -  w')^£  c to' 


G 


(88) 


where  r  and  7P  are,  respectively,  the  octahedral  shear  stress  and  plastic 
octahedral  shear  strain, 


w  2  Zq  -  zg  ,  (89) 

and  zg  denotes  the  value  of  Zq  at  the  end  of  the  pure  hydrostatic  compression 
phase  of  the  test.  In  Eq.  (88),  r(w)  and  d7P/dw(w)  are  presumed  known  while 
p( w)  is  to  be  determined.  To  determine  r(w)  and  d7P/dw(w),  the  relation 
between  7P  and  w  during  shearing  must  be  known.  This  relation  was  found  by 
numerically  integrating  Eq.  (49),  which  can  be  rewritten  in  the  form: 


33 


SSS-R-88-9375 


dw  = 


2*6  . 
e  +  5 

Q _ 

2 

d7P 

I'o  re 

e2"S  +  3 

- 1 

2 

•  1 

1/2' 

* 

(90) 


Inasmuch  as  both  r  and  cR  are  known  as  functions  of  7P  from  Eqs.  (83)  and 
(84),  the  righthand  side  of  Eq.  (90)  can  be  expressed  solely  in  terms  of  7",  i.e., 
Eq.  (90)  can  be  written  in  the  form: 


1 


dw  =  f  (7P)  d7P 


(91) 


and  numerically  integrated  to  give  the  relation 


w  ■  w(7P)  (92) 

In  this  manner,  the  relation  between  w  and  7P,  shown  in  Figure  8,  was  obtained. 
The  relation  (92)  was  numerically  inverted  to  give 


7P  =  7P(w)  (93) 

Therefore,  by  using  Eq.  (93),  the  relations  (83)  and  (84)  were  expressed  in  the 
form: 


r  =  r(w) 

(94) 

4?  -  $£<»> 

Using  these  relations  in  conjunction  with  Eq.  (88),  the  resulting  Volterra  integral 
equation  was  solved  numerically  for  i(w),  using  a  technique  described  in  Ref.  3. 
The  result  was  then  fit  by  Prony’s  method  [14]  to  yield  the  following  expression 
for  the  shear  kernel  function: 
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CM  -  /  V  *r2°  I®51 

r=1 

where 

k  =  11.6  N/Pa 
A 2  =  160  M’a 
Ag  -  706  IVPa 

To  provide  a  check  on  the  internal  consistency  of  the  above  operations,  we 
performed  the  following  analysis.  The  expression  given  by  Eq.  (95)  for  p(zk), 
together  with  the  derivative  d7wdw  obtained  from  the  relation  between  7P  ana  w 
depicted  in  Figure  8,  were  used  in  conjunction  with  Eq.  (88)  to  predict  r(w). 
Figure  9  shows  a  comparison  between  the  predicted  r(w)  and  the  coresponding 
r(w)  from  the  data  for  test  RDC  624.  As  the  figure  reveals,  the  agreement  is 
quite  good  and  thus  validates  the  accuracy  of  the  above  technique. 

To  determine  the  shear  hardening  function,  FD,  we  first  recall  that  in  the 
preceding  devlopments  we  set  Fq  =  1  at  the  reference  pressure  of 
<rR  =  6.90  MPa.  In  view  of  this,  we  can  write  Eq.  (70)  as 


a 1  *  9.0 

a2  =  28.9  (96) 

«3  *1,381. 


Fd(jo)  = 


(97) 


since  sw  and  rm  differ  only  by  a  multiplicative  constant.  The  dependence  of  rm 
on  aQ,  for  the  ISST  soil,  is  given  by  the  data  in  Figure  10  and  is  found  to  have  the 
following  linear  form: 


T.(ffo)  “  Ts  +  Vo  (98) 
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r  versus  7P  for  test  RDC  624,  showing  comparison  between 
model  prediction  and  data. 
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Figure  10.  rm  versus  a  for  ISST  soils.  The  circles  denote  data  points  from 
Ref.  4. 
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where 

rg  «  0.15  M^a 
Ps  -  0.56 

For  the  reference  pressure  of  6.9  MPa,  we  find  from  Figure  7  that 
r.(*o)  *  4.0  NPa 

Thus,  from  Eqs.  (97)  to  (100),  it  follows  that 
fD  ■  'o  +  V 

where 

*0  =  0.0375 
^  =  0.140  VPa'1 

The  specification  of  the  model  parameters  for  the  ISST  soil 
complete. 
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Section  4 

NUMERICAL  CONSIDERATIONS 

A  numerical  scheme  is  presented  in  this  section  for  integrating  the  system 
of  equations  which  govern  the  new  endochronic  model  with  jpilatancy  described 
in  Section  2.2.  The  scheme  is  explicit,  first  order  accurate  and  based  upon 
Euler's  method.  Accordingly,  care  must  be  taken  in  applying  the  method  to 
ensure  that  the  prescribed  increments  are  sufficiently  small  so  that  the  computed 
results  are  essentially  independent  of  the  increment  size.  Otherwise,  the  scheme 
is  straightforward,  efficient  and  easy  to  implement. 

Two  different  versions  of  the  scheme  are  described  below,  namely,  one 
which  requires  the  prescribed  strain  history  as  input  and  one  which  requires  the 
prescribed  stress  history  as  input.  In  both  cases,  the  equations  are  restricted  to 
principal  directions  of  stress  and  strain.  Another  version  that  specifically  applies 
to  the  standard  triaxial  compression  configuration  has  been  developed  under 
separate  contract  [16].  A  listing  of  a  computer  program  for  the  version  of  the 
scheme  that  requires  the  prescribed  strain  history  as  input  is  given  in 
Appendix  B. 

4.1  BASIC  EQUATIONS. 

The  basic  equetions  which  govern  the  new  endochronic  soil  model  with 
dilatancy  were  giver*  earlier  in  Section  2.2  but  are  summarized  below  for 
convenience; 


§  = 


rft> 


os” 

z’)ar 


dr 


(103) 


a  = 


rpH  ■  z')ar dz'  + 


rpH  - z')  s  *  ar dz' 


(104) 


A  second-order  accurate  scheme  has  recently  been  developed  by 
Murakami  and  Read  [15]. 
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i 


dg  =  ae[d9  -  dgP) 

(105) 

da  =  «(de  -  dep) 

(106) 

dz2=  ||dgP||2  +  k2(d<P)2 

(107) 

dip  =  dz/Fp 

(108) 

d^  -  dz/  (kF^ 

(109) 

Here,  g  denotes  the  deviatoric  stress  tensor,  a  is  the  hydrostatic  stress 
(pressure),  gp  represents  the  plastic  component  of  the  deviatoric  strain  tensor  g, 
while  ep  is  the  plastic  component  of  the  volumetric  strain  e.  Moreover,  G  and  K 
are,  respectively,  the  shear  and  bulk  moduli,  while  k  is  a  contant  which 
determines  the  magnitude  of  shear-volumetric  coupling.  The  double  bars 
surrounding  a  symbol  denote  its  Euclidean  norm,  while  single  bars  denote 
absolute  value.  Furthermore,  Fq  and  Fu  are,  respectively,  shear  and  hydrostatic 
hardening  functions.  In  addition,  z  denotes  the  intrinsic  time  scale,  while  Zq  and 
Zj_l  are,  respectively,  the  intrinsic  times  for  shear,  dilatancy  and  hydrostatic 
behavior.  Finally,  p(z)t  ^(z)  and  T(z)  are  weakly  singular  kernel  functions 
satisfying  the  condition  /»(0)  =  ^(0)  =  T(0)  =  »,  but  integrable  in  the  domain 
0  i  z  <  •. 

The  weakly  singular  kernel  functions  can  be  expressed  in  terms  of  Dirichlet 
series: 


I 

I 


P(z) 


i 

*(Z) 


(110) 


(HI) 


i 


l 

l 

i 


! 


I 
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(112) 


whore  in  order  to  satisfy  the  Clausius-Duhem  inequality,  it  is  necessary  that 
ar  1 0, *  0, 7j  I  0  and  i  0,  2  0,  l*j  2  0.  Moreover,  to  ensure  that  f(z),  f(z) 

and  r(z)  are  singular  at  tne  origin  and  integrable  over  a  finite  domain,  we  must 
have 


4.2  INCREMENTAL  FORM  OF  BASIC  EQUATIONS. 

In  applications  of  the  theory  to  date,  it  has  been  found  that  two  or  three 
terms  of  the  series  (110)  to  (112)  are  usually  quite  adequate  for  representing  the 
kernel  functions.  In  such  cases,  however,  care  must  be  taken  to  ensure  that  the 
infinitely  large  values  of  p( 0),  1(0)  and  r(0)  are  approximated  by  suitably  large 
finite  values.  When  this  is  done,  we  can  write: 


(115) 
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(116) 


where  n  and  m  are  finite. 

It  then  follows  that  the  expressions  for  g  and  e,  as  given  by  Eqs.  (103)  and 
(104),  can  be  alternately  written  as 


(117) 


(118) 


where  Qj,  P,  and  Nj  satisfy  the  following  ordinary  differential  equations: 


From  Eqs.  (117)  to  (121)  we  can  write: 


(119) 

(120) 

(121) 
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dj-AdjP  -9*0 

(122) 

d»  -  B  d«P  +  r(*  •  c^P) 

•  (P  +  N)dz^ 

(123) 

where 

m 

N  -  N, 

(124) 

m 

P  -  P, 

Equations  (122)  and  (123)  provide  a  simple  approach  for  incrementally 
integrating  the  stresses  $  and  cr,  which  is  considerably  more  attractive  from  a 
computational  standpoint  than  numerically  coping  with  the  hereditary  integrals  in 
Eqs.  (103)  and  (104). 

In  that  which  follows,  explicit  numerical  schemes  are  presented  for 
incrementally  updating  the  endochronic  equations  given  above  when  either  the 
strain  incrments  or  the  stress  increments  are  given.  Because  of  the  explicit 
nature  of  the  scheme,  it  is  necessary  that  the  increments  be  taken  sufficiently 
small  to  ensure  accuracy. 

4.3  PRESCRIBED  STRAIN  INCREMENTS  A*. 

It  is  assumed  that  g,  o,  g,  gP,  e,  eP,  Q.,  Nj  and  Pj  are  known  at  the 
beginning  of  each  prescribed  strain  incrment,  Aj.  From  Eqs.  (108),  (109),  (122) 
and  (123),  we  can  write 
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4*-AA»P-0f* 

(125) 

to  *  B  AeP  +  r(s  •  AgP)  -  P"[y-^1  Az 

I  H  J 

(126) 

If  we  now  combine  these  equations  with  the  incremental  Hooke's  Law, 
Eqs.  (105)  and  (106),  it  follows  that 


*8P  -  -pr^W  (®  *» +  o  IQ  <127> 


4tP- 


TETF-RT 


(K&6  - 


Cffe  (s 


*8) 


P  +  N 
H 


r(0*s?  y 

CT5|AzJ 


(128) 

Upon  substituting  these  results  into  Eq.  (107),  the  following  quadratic  expression 
for  Az  is  obtained: 


a  Az*  +  b  Az  +  c  -  0 


(129) 


where 


a  =  1  - 


9  •  9 


(A+2S)2f£  (B  +  K)2 


P  +  N  ^(9  *  j) 


"H  (A  +  23)  Fq 


(130) 


b  =  -  2 


f  2GQ*  AS 


[(A +23)%  (B  +  K)‘ 


P  +  N  rIP 

'  T^rssTFp 


^S5  (S  *  ‘8)  ■  K&e 


(131) 
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<=  ■  •  {br?y * 4*  ’  +  ^^3  K  •  jftfc  (» •  4*ll2) 

(132) 


it  can  be  verified  that  the  above  expression  for  a,  b,  and  c  reduce  to  those  given 
In  Ref.  17  when  dilatancy  is  not  present,  i.e.,  N  »  I*  -  0.  Equation  (129) 
provides  two  roots  Az1  9,  the  one  of  interest  being  the  one  for  which  Az  *  0. 
Once  Az  is  known,  Ag^  And  A«^  can  be  obtained  from  Eqs.  (127)  and  (126),  after 
which  A$  and  A*  can  be  obtained  from  Eqs.  (125)  and  (126).  Finally,  Eqs.  (119) 
to  (121)  are  used  to  update  the  Nj  and  P|.  This  approach,  therefore,  permits 
one  to  determine  the  stress  increments,  A£,  for  prescribed  increments  in  the 
strain  Ag. 

4.4  PRESCRIBED  STRESS  INCREMENTS  l(. 

In  this  case,  5. »,  g,  jp,  e.  «p.  QL  N,  and  P,  are  assumed  to  be  known  at  the 
beginning  of  each  prescribed  stress  Increment  Ag.  From  Eqs.  (109),  (122)  and 
(123).  we  can  obtain  the  expressions: 


A#" 


As  +  p-  A* 
hD 


“P  -  B 


*»  •  £  (S  *  *S)  + 


p  +  N  r  (9  *  *) 

.■wsr '  * 


Az 


(133) 


(134) 


Upon  substituting  these  results  into  Eq.  (107),  we  obtain  the  following  quadratic 
expression  for  Az: 


a  Az2  +  b  Az  +  c  «  0 
where 


(135) 


46 


SSS-R-88-9375 


a  =  1 


Q’Q 

(AFo)2 


1_ 

B2 


8] 


X Tv 


(136) 


b  = 


-  a  [s  *  ‘si 


P  +  N  _  r(9  *  §)1 
^FH  A  fd 


(137) 


(138) 


Again,  the  root  of  interest  from  (*q.  (135)  is  that  for  which  Az  l  0.  Once  Az  is 
found  from  Eq.  (135),  Ag  and  Aep  can  be  determined  from  (133)  and  (134),  after 
which  Ag  and  Ac  may  be  found  from  Eqs.  (105)  and  (106).  The  Qr,  Nj  and  the  Pj 
are  updated  on  the  basis  of  Eqs.  (119)  to  (121).  This  approach,  therefore, 
permits  one  to  determine  the  strain  increments  Ag  for  prescribed  increments  in 
the  stress  Ag. 
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Section  5 


NUMERICAL  ANALYSIS  OF  SHEAR  AT  FIXED  HYDROSTATIC  PRESSURE 

In  order  to  provide  a  check  on  the  accuracy  of  (1)  the  numerical  scheme, 
(2)  the  computer  program,  and  (3)  the  material  functions  and  parameters  that 
had  been  determined  for  ISST  soil,  we  conducted  a  computer  study  of  the  four 
test  cases  of  shear  at  fixed  hydrostatic  pressure  reported  by  WES  in 
Reference  A.  The  stress-drive  version  of  the  numerical  scheme,  described  in 
Section  4.4,  was  used  for  this  purpose,  since  the  corresponding  laboratory  tests 
had  been  conducted  under  stress-controlled  conditions.  The  general  stress 
path  followed  by  these  tests  in  the  Rendulic  plane  is  depicted  in  Figure  1 1 ,  which 
shows  the  four  legs  consisting  of  (1)  pure  hydrostatic  loading  up  to  some  fixed 
hydrostatic  pressure  aQ,  (2)  shear  loading  at  fixed  aQ,  (3)  shear  unloading  at 
fixed  a Q,  and  (4)  pure  hydrostatic  unloading. 

5.1  DESCRIPTION  OF  DIFFICULTY. 

The  numerical  calculations  proceeded  smoothly  and  gave  satisfactory 
results  for  the  first  three  legs.  However,  on  the  fourth  leg,  which  involved  pure 
hydrostatic  unloading,  difficulties  typically  arose  after  the  hydrostatic  pressure 
had  been  reduced  from  between  20  to  25  percent  of  its  peak  value  aQ.  To 
elaborate  on  the  difficulty  in  greater  detail,  let  us  first  recall  that,  in  the  case  of  the 
stress-driven  numerical  algorithm,  the  increment  in  the  intrinsic  time  Az  is 
determined  from  the  quadratic  expression  given  by  Eq.  (135),  i.e., 


a  Az2  +  b  Az  +  c  =  0  (139) 

where  a,  b  and  c  are  defined  by  Eqs.  (136)  to  (138),  but  may  be  alternately 
expressed  in  the  following  form: 


(140) 
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Description  of  stress  path  in  Rendulic  plane  for  WES  tests  of  shear 
at  fixed  hydrostatic  pressure. 


49 


SSS-R-88-9375 


k  P* 

ZTh 


(141) 


(142) 


if  we  set 


La*  =  Aa  -  t-  (§  •  A§) 
H 


(143) 


P*=P-^(S*0) 


(144) 


Inasmuch  as  Az  is,  by  definition,  unique  and  positive,  an  admissible  solution  of 
Eq.  (139)  is  obtained  when  the  two  roots,  Az,  and  Az2,  are  real,  unequal  and  of 
opposite  sign.  When  this  is  the  case,  Az  is  defined  by  the  positive  root,  and  the 
negative  root  is  considered  redundant.  An  admissible  solution  is  therefore 
obtained  when  the  following  inequalities  are  satisfied: 


b  -  4ac  >  0 


(145) 


c/a  <  0 


(146) 


Note  from  Eq.  (142)  that  the  parameter  c  is  always  negative.  The  parameter  b 
may  be  either  positive  or  negative.  Thus,  if  a  <  0  then  either  real  roots  Az  do  not 
exist,  or  they  both  have  the  same  sign.  Either  of  these  two  cases  is  inadmissible 
and  the  computer  program  was  designed  to  stop  if  this  situation  arose.  The 
parameter  a  must  be  positive  if  an  admissible  solution  is  to  be  obtained. 
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For  the  cases  of  shear  at  constant  hydrostatic  pressure  that  we  studied 
numerically,  the  parameter  a  became  negative  in  the  early  stages  of  the 
hydrostatic  unloading  leg  and,  when  this  occurred,  the  computer  program 
stopped.  To  obtain  an  understanding  of  how  this  can  occur,  let  us  examine  the 
governing  equation  for  Az  on  this  leg.  Note  that  for  hydrostatic  unloading  under 
zero  shear,  A§  =  0  so  that  A  a*  =  A  a.  In  this  case,  the  solution  of  the  quadratic 
equation  (139)  yields  the  following  result: 


Az 


k  A  a 

p*  „ 

1  - 

IT’ 

"2 

B 

i  pi 

1  '  51 

.  » 

a  ■ 

IT 

2 

• 

(147) 


while  the  parameter  a  in  Eq.  (140)  becomes 


If  Fq  goes  to  zero  as  a  tends  to  zero,  as  is  the  case  with  ISST  soil  if  we  disregard 
a  very  small  cohesive  strength,  then  it  is  possible  for  a  to  become  negative,  i.e., 


1  - 


0 

A  fd 


P*  12 

TO 


<  0 


or 


(149) 


1  - 


9 


P*  i; 


(150) 


Thus,  in  view  of  Eq.  (147),  both  roots  of  Az  will  have  the  same  sign. 
Furthermore,  since  A  a  is  negative  during  hydrostatic  unloading,  both  roots  will  be 
positive  and  therefore  inadmissible  since  one  does  not  know  which  of  the  two 
positive  roots  is  the  correct  one,  i.e.,  there  is  a  uniqueness  problem. 
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5.2  SOURCE  OF  DIFFICULTY. 

In  order  to  determine  whether  or  not  the  difficulty  is  of  numerical  origin,  we 
took  the  following  steps: 

1.  First,  each  of  the  problems  was  rerun  with  a  variety  of  different  size  stress 
increments,  which  varied  by  over  an  order  of  magnitude.  As  an  extreme 
case,  we  used  20,000  stress  increments  on  the  fourth  leg  alone. 
Regardless  of  the  size  of  the  stress  increment,  the  computer  program 
stopped  (because  of  the  condition  a  <  0)  at  virtually  the  same  identical 
value  of  the  hydrostatic  pressure  on  the  unloading  leg  in  each  case.  The 
magnitude  of  the  hydrostatic  pressure  at  which  this  occurred  was  different, 
however,  for  each  of  the  four  cases  considered.  Thus,  the  difficulty  did  not 
appear  to  be  due  to  increment  size. 

2.  Secondly,  we  replaced  the  usual  approach  for  obtaining  the  roots  &z  of  the 
quadratic  equation  (139)  in  the  computer  program  by  a  more  accurate 
approach.  As  most  numericists  are  aware,  if  the  usual  expressions  are 
used  to  determine  the  roots  of  a  quadratic,  serious  difficulties  can  arise.  In 
particular,  if  either  a  or  c  (or  both)  are  small,  then  one  of  the  roots  will 
involve  the  subtraction  of  b  from  a  very  nearly  equal  quantity  (the 
discriminant)  and  the  root  will  be  determined  very  inaccurately.  The  correct 
way  to  compute  the  roots  is  as  follows.  If  we  set 


q  = 


b  +  sgn(b)Jb2  -  4  ac 


(151) 


then  the  two  roots  are 


Az1  =  9  and  Az2  =  §  (152) 

This  method  of  determining  the  roots  &z1  and  Az„  was  introduced  into  the 
computer  program  and  had  no  perceptible  effect  on  the  result*.  The 
computer  calculations  continued  to  stop  at  the  same  locations  on  the 
hydrostatic  unloading  legs. 
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Next,  the  method  of  integrating  the  differential  equations  (119)  to  (121)  for 
the  Qr,  Pj  and  was  improved.  Consider,  for  example,  the  differential 
equation  for  Q  given  by  Eq.  (119).  A  careful  analysis  shows  that  this 
differential  equation  can  be  integrated  over  a  finite  increment  Az  to  yield  the 
following  expression  for  the  inu  ament  AQr: 


(153) 


where  AZq  *  Az/Fg.  Similar  expressions  are  obtained  for  AP,  and  AN-. 
Inasmuch  as  the  term  inside  the  first  parentheses  on  the  right  side  of  this 
equation  has  the  potential  for  numerical  difficulties  when  Azp  is  small,  we 
replaced  this  term  with  an  accurate  Pade  approximations  avoid  the 
difficulty.  Again,  the  introduction  of  this  improvement  into  the  computer 
program  had  no  perceptible  effect  on  the  calculational  results  and  the 
calculations  continued  to  stop  at  nearly  the  same  location,  as  previously, 
on  the  hydrostatic  unloading  leg. 


In  view  of  the  fact  that  none  of  the  changes  described  above  to  the 
numerical  scheme  had  any  significant  effect  in  alleviating  the  problem  of  the 
parameter  a  becoming  negative  on  the  hydrostatic  unloading  leg,  it  was 
concluded  that  the  problem  was  not  numerical  in  nature.  Attention  was  then 
turned  to  the  model  itself. 


We  initially  speculated  that  the  problem  may  arise  from  the  Dirichlet  series 
representations  for  ^(Zq)  and  f(zH)  anc|.  t0  ®xP,ore  this.  we  developed 
approximate,  but  reasonable,  representations  of  ^(zD)  and  ^(zH)  based  on  the 
use  of  a  single  exponential  in  each  case,  which  corresponds  to  one  internal 
variable.  Using  these  representations  for  p{ zD)  and  #(zH),  numerical  studies  of 
the  problems  of  shear  at  constant  pressure  were  redone,  but  this  time  the 
calculations  proceeded  to  completion  without  difficulty  over  the  entire  four  legs 
of  the  stress  path.  These  results  would  certainly  indicate  that  the  difficulty  was 
due  to  the  Dirichlet  series  representations  for  p(zp)  and  #(Zu).  On  further 
reflection,  however,  it  was  concluded  that  this  is  not  realistic.  Consider  that  a 
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Oirichlet  series  represents  a  solid  of  some  kind  and  the  only  thermodynamic 
constraint  is  that  the  series  consist  of  positive  decaying  exponential  terms.  If  we 
accepted  the  above  plausibility,  one  must  conclude  that  some  solids  suffer  from 
this  difficulty  and  some  do  not,  a  conclusion  which  is  physically  unrealistic. 

We  now  believe  that  the  problem  lies  in  the  physics  of  the  representation  of 
the  material  behavior.  For  instance,  in  the  developments  to  date,  the  parameters 
Aj.  have  been  assumed  constant  and  the  hardening  function  Fp  has  been  taken 
to  be  a  function  of  o  only,  although  there  is  strong  evidence  that  it  also  depends 
on  e^.  To  elaborate,  we  note  that  various  soils  with  different  initial  porosities 
exhibit  different  responses  to  shear.  Hence,  porosity  must  affect  Fp.  However, 
since  porosity  and  plastic  volumetric  strain  are  interrelated,  plastic  volumetric 
strain  must  also  affect  Fq.  It  is  very  likely  that  the  same  will  be  true  of  the 
parameters  A  .  A  study  is  needed  to  examine  this  dependence  with  a  view 
toward  ensuring  that  the  parameter  "a"  in  Eq.  (140)  is  never  negative,  whatever 
the  stress  (or  strain)  path.  This  will  ensure  that  the  difficulty  will  not  arise. 


Section  6 


CONCLUSIONS  AND  RECOMMENDATIONS 

On  the  basis  of  the  analysis  described  in  preceding  sections  of  this 
report,  the  following  conclusions  and  recommendations  are  given: 

1.  A  new  endochronic  plasticity  model  for  soils,  which  can 
describe  both  densification  and  dilatancv.  has  been 
successfully  applied  to  laboratory  data  generated  by  WES  for 
ISST  soils. 

2.  Appropriate  theoretical  methods  have  been  developed  from 
which  the  material  functions  and  parameters  of  the  model  can 
be  evaluated  directly  from  two  types  of  soil  tests,  namely, 

(a)  pure  hydrostatic  compression  tests  and  (b)  triaxial  tests  in 
which  the  soil  is  sheared  at  fixed  hydrostatic  (not  confining) 
pressure. 

3.  An  explicit,  efficient  numerical  scheme  was  developed  for 
numerically  integrating  the  system  of  equations  which  govern 
the  model.  A  corresponding  computer  program  based  on  the 
numerical  scheme  was  also  developed  and  applied  to  several 
problems. 

4.  A  difficulty  with  the  governing  equations  arose  when  the 
computer  program  was  applied  to  the  case  of  shear  at 
constant  hydrostatic  pressure.  Specifically,  on  the  final  leg  of 
the  stress  path  for  this  case,  which  involves  pure  hydrostatic 
unloading,  the  computer  program  calculated  inadmissible 
increments  in  the  intrinsic  time  increment  Az. 

5.  A  number  of  aspects  of  the  numerical  scheme  were  improved 
but  none  appeared  to  have  any  significant  effect  on  alleviating 
the  difficulty  noted  above.  Decreasing  the  increment  size  by 
over  an  order  of  magnitude  also  had  no  noticeable  effect.  In 
view  of  this,  it  was  concluded  that  the  difficulty  was  not  of 
numerical  origin. 
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6.  We  believe  that  the  difficulty  lies  in  the  physics  of  the 
representation  of  the  material  behavior.  It  appears  that  the 
mathematical  representations  adopted  in  the  present  model 
for  some  of  the  material  functions  may  not  be  general  enough 
to  describe  some  of  the  observed  behavior  of  soils.  A  study 
is  needed  to  explore  this  issue,  with  the  goal  of  ensuring  that 
the  difficulty  noted  above  will  not  arise,  whatever  the  stress  or 
strain  path. 

In  closing,  we  emphasize  that  this  study  represents  the  first  attempt  to 
apply  the  new  endochronic  plasticity  model  with  dilatant  capability  to  real  soils 
and,  in  the  course  of  the  study,  considerable  new  insight  into  the  characteristics 
of  the  model  was  obtained.  In  view  of  the  substantial  promise  shown  by  the 
model  in  describing  the  complex  features  of  real  soil  behavior,  we  strongly 
recommend  that  a  further  study  be  undertaken  to  resolve  the  issue  of  the 
representation  of  material  functions  discussed  above. 
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Appendix  A 

TRANSFORMATION  FROM  ENGINEERING 
STRAIN  TO  NATURAL  STRAIN 


The  ISST  soil  data  racantly  supplied  by  WES  to  S-CUBED  for  modal 
development  are  expressed  in  terms  of  engineering  strains.  Most  of  the 
materials  models  used  in  the  defense  community,  however,  are  based  upon  the 
natural  (or  logarithmic)  definition  of  strain.  In  developing  a  constitutive  model 
from  the  WES  data  for  eventual  use  in  conjunction  with  one  of  the  material 
response  coeds,  the  engineering  strains  given  in  the  data  must  be  transformed 
to  natural  strains  for  consistency.  The  purpose  of  this  appendix  is  to  document 
the  equations  for  making  this  transformation. 

Analysis. 

Consider  the  small  element  of  material,  shown  in  Figure  1.  which  has  sides 
of  length  L°,  l£,  and  in  the  initial  (unstrained)  state. 


zz 

7 

A 

c 


Figure  1.  Unstrained  configuration  of  material  element 


The  element  is  strained  in  (principal)  directions  perpendicular  to  its  surfaces  so 
that  the  lengths  of  thesides  become  L1 ,  Lg  and  Lg.  Under  these  conditions,  the 
engineering  strains,  e  j ,  can  be  expressed  as 


(i  =  1,2,3) 


while  the  corresponding  natural  strains,  e  j ,  are  given  by 


(D 
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The  minus  signs  appear  on  the  right-hand  sides  of  Eqs.  (i)  and  (2)  so  that  the 
strains  wit  be  positive  in  compression,  which  is  the  sign  convention  commonly 
adopted  in  sofi  mechanics. 

Equation  (1)  can  be  solved  for  Lj  to  give 


Li”L?(1**Ll  ('  *1.2.3)  O) 

which,  when  substituted  into  Eq.  2,  leads  to  the  expression 


This  equation,  therefore,  relates  the  engineering  strains  to  the  natural  strains. 

Let  us  consider  now  the  volumetric  strain.  The  engineering  definition  of 
volumetric  strain  is 


V  -  V. 


(5) 


where  V  and  vQ  denote,  respectively,  the  current  and  initial  volumes  of  the 
element  shown  in  Figure  1 .  In  terms  of  natural  strain,  the  volumetric  strain  is 
given  by 
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(6) 


In  terms  of  the  current  and  initial  lengths  of  the  sides  of  the  small  element,  we 
can  write 


V 

VI  i  o.  o,  o 

Equation  (3)  may  be  used  in  Eq.  (7)  to  give 


Combining  this  with  Eqs.  (5)  and  (6),  we  find 


(7) 


(8) 


and 
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b  -  4) 

] 

(9) 


(10) 


By  setting 
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'e  =  ‘?+«|+‘3 

"c  ■  «l'i  +  ll*l  +  *3*1  (”) 

,Mt“el'2e3  • 

Equations  (9)  and  (10)  can  be  re-expressed  in  the  forms: 

*v~  •«  -  "e  ^ IM*  <«) 

*v“  '  ***(1  ''  I*  +  "*  *  <13> 

For  small  strains,  Eqs.  (12)  and  (13)  reduce  to  the  usual  expressions  for  the 
volumetric  strain  given  by  the  infinitesimal  theory,  namely, 


(14) 


Summary. 

Equation  (4)  provides  the  transformation  from  axial  engineering  strains 
to  natural  strains  e  j ,  while  Eq.  (10)  allows  ope  to  obtain  the  natural  volumetric 
strain  from  the  axial  engineering  strains,  e j . 
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Appendix  B 


COMPUTER  PROGRAM  FOR 
ENDOCHRONIC  SOIL  MODEL, 
REQUIRING  THE  STRAIN  HISTORY 
AS  INPUT. 
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Index  convention: 
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*  *  l  3  l 

,  sig  (  3  ) 

(maxi) 

,  biggamfmaxi) 

,  beta(maxi)  , 

(maxi) 

,  am  1  gam  (maxi) 

9 

(maxr) 

,  alpha  (maxr) 

,  qr(3,maxr)  , 

(maxi) 

,  epsleg(3,max 

1) 

ds  ( 
epso  ( 


1: 


bign(maxi) 


read 

read 

read 

read 

read 

read 

read 


210, 

210, 

*. 


nof i I e 
ndf i  le 

akzz  ,am  ,amuz  , both , bets , betk 
capkz,capkl, const, pr  ,taus,cO  ,sigt 
numr 


numi 


r  ,  (a(j)  ,alpha(j)  ,  j=l,numr) 

i  ,(b(i),beta  (i) ,biggam(i) ,smlgam(i) , i=l,numi) 


*,  numleg,nleg(l) , (epsleg(k,2) ,k=l,3) 


print  220,  akzz  ,am  ,amuz  , both, bets, betk 
print  230,  capkz,capkl, const, pr  ,taus,cO  ,sigt 
print  240,  (j,a  (j),j,»lpha  (j) , j=l,numr) 

print  250,  (i,b  (i),i,beta  ( i) , i=l,numi) 

print  260,  (i ,biggam(i) , i ,smlgam(i) , i=l,numi) 
print  270,  (l,nleg(l) , (epsleg(k,l+l) ,ksl,3) , 1=1,1) 


open  f uni  tall  ,fi  l< 
open (uni t=12,fi l< 
read  (11,*) 
read  (11,*) 
read  (11 ,*) 
road  (Hi*) 
read  (11,*) 
read  (11,*) 
read  (11,*) 
read  (11,*) 
read  (11,*) 

write  (12,310) 

capasO. 

do  10  j&l,numr 
10  capaacapa+a(j) 

capbaO. 

gamOaO. 

do  20  i=l,numi 


endf i I e, forms* formatted*) 
snof i le,forma*formatted*) 
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capb*capb*b(i) 

20  gaaO»flaiB(kb  i  ggaai  ( i ) 

n«l 

do  120  nlal,maal«g 

opoldlaopnowl 
epo I d2*opnow2 
opo I d3*opnow3 

road  (11,*)  opnowl , opnow2 , opnowS 

opnowla-epnowl 

opnow2w-opnow2 

opnow3»-opnow3 

do  110  nnsljnlog(l) 

nstopssn 
n  «n+l 
opso(l)*ops(l) 

•pso  (2)  sops  (2) 
epsoi3)seps(3) 

ops  ill sepo I dl+ (opnowl -opo I dl) * (f I  oat  f nn) ) /f I  oat (n I og  (1) ) 
ops  i2 1 atopo  I  d2+  (epnew2-opo  I  d2 1  *  if  I  oat  f  nn  1 )  /f  I  oat  (n  I  og  (1)  ) 
ops  (3)  «opo  I  d3+  (opnew3-opo  I  d3)  *  (f  I  oat  (nn)  )  /f  I  oat  (n  I  og  (l)  ) 

dov«0. 

do  40  k*l,3 

dops (k) sops (k) -opso (k) 

40  dovadov-tdops(k) 

•v  sov  +dov 

do  60  k*l,3 
q(k)-0. 

do  50  j*l,numr 

50  q  (k)*q  (kUalpha(j)*qr(k, j) 
o  (k)aops  (k)  —  ov/3. 

60  do(k)sdeps(k)-dov/3. 

p  =0. 

f  n=0. 

do  70  isl^umi 
f  n=fn+smlgam(i)*bign(i) 

70  p  =p+bota(i)*pi (i) 

qde  =0. 
qq  =0. 
qs  =0 . 
sde  =0. 
dodesO. 
do  80  k*l,3 
qdo  *qdo  ♦q  (k)*de(k) 
qq  »qq  *q  (k)*q  fk) 
qs  sqs  +s  (ki*q  (k) 
sdo  *sdo  +s  (k)*do(ki 
80  dodo*dodo+do(k)*do(k) 
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Jfc-rt  &  *  K  A  An  1 .1  lAiBLTl  M'S  u 


e 


c 


c 


c 


c 


akz  aakzz*dexp(-betk*evp) 

capkacapkz*capkl*pp 

aau2«aau»Kow$t«d»qrt  (gaaz) 

gaa  agaaO*exp(-cO*evp) 

fa  ■( taus«b«ta*pp)  /  (taua«bata*pr) 

fh  adexp(beth*evp) 


abaracapa*a»u2 

bbaracapb+capk 

am  I aal . -qq/ (abar*f a) **2 

- (1 . /bbar**2) • ( (p+f  n) /f h-akz*gan*qa/ (abar*f a) ) **2 
am  I ba-2* (aau2*qde/ (abar **2*f a) ♦ (akz/bba  iO  **2 

* (capk*dev-gae*aeu2/abar*ade) * ( (p+f n) / (akz*f h) 

-gam*qa/ ( abar *f a) ) ) 

am I c«- ( (amu2/aba  r) **2*dede+ (akz/bba r) * *2* (capk*dev 
-  anu2*gaa*ade/abar) **2) 
arg  aaarib**2-4.*aala*amlc 

dzl  ■-(.5/smla)*(anlb<*>dsign(l.d0lsmlb)*daqrt(max(0.targ))) 
dz2  *smlc/(smla*dzl) 
if  (arg. It.O. .or .dzl*dz2.gt.0.)  then 
print  290,  nt ,natepa,arg,dzl,dz2 
end  if 

dz  *max(dzl,dz2) 

dzsadz/fs 

dzhadz/(akz*fh) 

devp»(l . /bbar) * (capk*dev-anu2*gan 

/abar*ade+ ( (p+f n) / (akz*fh) -g»a*q*/ (aba  r*f a) ) *dz) 
evp  aevp+devp 
dwp  mO. 
do  90  k*l,3 

dep  (k)  «  (amu2*de  (k)  *q  (k)  *dz/f a)  /aba  r 

ep  (kiaep(k)  +dep(k) 

da  (k) acapa*dep (k) -q (k) *dz/f a 

a  (k;  as  (k)+da(k) 

dwp  adwp+a (k) *dep (k) 

do  90  jsl,nunr 

emaladexp  f-alpha(j)*dzs) 

ems2adxpldx(  alpha(j)*dzs) 

ems  aemal*ema2 

90  qr (k , j ) =qr (k , j ) +ema* (a ( j ) *dep (k) -a  I pha ( j ) *qr (k , j ) *dzs) 
ppspp+capb*devp+gam*dwp- ( (p+f n) / (akz*f h) ) *dz 
do  925  k=l,3 
925  aig(k)ss(k)+pp 

do  100  ial,numi 

b i gn ( i ) =b i gn ( i ) +gamO*exp (-cO*evp) *dwp 
-b i gn ( i ) *a« I  gam ( i ) *dz/ (akz*f h) 
emhadexp (-beta ( i ) *dzh) adxpldx (beta ( i ) *dzh) 

100  pi  (i)spi  (i)+emh*fb(i)*devp-beta(i)*pi  (i)*dzh) 
gamz  adaqrt(fepa(l)->epa(2))**2 

♦  (epa(2)-eps(3) )**2 

♦  (epa(3;-epa(lj)**2)/3. 

if  (aig(l) . le.sigt.or.aig(2) . le.sigt.or .aig(3) . le.sigt)  go  to  130 
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writ#  (12,320)  n,(eps(k),kal,3), (aig(k)  ,k»l,3) 
110  continue 
120  continue 
130  cal  I  exit 


210  format (a) 
220  format( 


230  format ( 


240  format (/, 
250  format (/, 
260  format 


akzz 

amuz 

■*,lpel3.6,  * 

-»,lpel3.6,/. 

am 

a*, IpelS. 6 

both 

betk 

«*,lpel3.6,  » 

-*,lpel3.6) 

bets 

a*, IpelS. 6 

capkz 

const 

■*,lpel3.6,  * 

-*,lpelS.6,/, 

capkl 

a* t IpelS. 6 

pr 

-»,lpel3.6,  * 

tauz 

a*,lpel3.6 

cO 

a*, IpelS. 6,  * 

sigt 

■*,lpel3.6; 

a»,il. 

*  ■*, lpel3.fi, * 

alpha*, il,*  a* 

smlgam*,il, 'a*, 
sigleg  **,lp3el^ 


»  b’,il,’  a»,lpel3.6,* 

_ _ _  v*  biggam*,  il, *a*,lpel3.6,  * 

270  format!/, (*  leg*, i2, *, *, i6, *  points, 

280  format! *nateps  **,i5) 

290  format! *nl ,nsteps,arg,dzl,dz2’,2i6,lp3el5.7) 

310  format!*  strain  drive  endochronic  model  results*,/, 

el  e2  e3 

*  sigl  sig2  sig3  *) 

320  format(i5,3fl0.6,3fl0.4) 


IpelS . 6) ) 
1 pel3. 6)) 
lpel3.6)) 
•6)) 
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